

************************************
*Numeric Export disabled
local TEX ""
local TEXON 0 //Change "0" to "1" to enable exporting numeric results to .txt files
if !`TEXON' local TEX "*"
************************************
************************************
*Figure Export disabled
local FIG ""
local FIGON 0 //Change "0" to "1" to enable exporting figures (RELEVANT FOLDERS MUST FIRST BE CREATED)
if !`FIGON' local FIG "*"
************************************

**************************************************************************************
**************************************************************************************
*NEED TO SET PATH TO WHERE "Enns_SociologicalScience_Reproduction" FOLDER IS SAVED
cd ""
**************************************************************************************
**************************************************************************************


*Load Gilens' data and recode variables
do "Enns_SociologicalScience_Reproduction\LoadGilensDataRecodeConfirmAccuracy.do"

**********************
*Correlations across preferences
cor pred10_sw pred50_sw pred90_sw
mat def cor = r(C)
local corlowhigh = cor[3,1]
local cormidhigh = cor[3,2]
cor pred50_sw PREDALL_SW
mat def cor = r(C)
local cormidpublic = cor[2,1]
******************************


************************************************
**Exact replication of Gilens 2012, Table 3.3 50th and 90th Income Percentiles, Row 3 (p.79)
************************************************
************************************************

************************************************
*Table 1, Table 2, Table A-1 (top half): Columns 1 and 3 (Gilens Replication)
************************************************
*************************
*10 vs 90 (ten vs ninety: tvn)
*************************
*10
logit policyadopted pred10_sw_lor  if  abs(pred90_sw - pred10_sw)>.1
*Store results as local vars
mat def results = r(table)
local coeflowtvn = results[1,1]
local selow = results[2,1]
local lowblowtvn = results[5,1]
local upblowtvn = results[6,1]
local coeflowtvncons = results[1,2]
local lowblowtvncons = results[5,2]
local upblowtvncons = results[6,2]
local ntvn = e(N)
****************
*90
logit policyadopted pred90_sw_lor  if  abs(pred90_sw - pred10_sw)>.1
mat def results = r(table)
local coefhightvn = results[1,1]
local sehightvn = results[2,1]
local lowbhightvn = results[5,1]
local upbhightvn = results[6,1]
local coefhightvncons = results[1,2]
local lowbhightvncons = results[5,2]
local upbhightvncons = results[6,2]

************************************************
*Table A-1 (bottom half), Table A-2: Columns 1 and 3 (Gilens Replication)
************************************************
*************************
*50th income percentile vs 90th income percentile
*************************
*50
logit policyadopted pred50_sw_lor if abs(pred90_sw - pred50_sw)>.1
*Store results as local vars
mat def results = r(table)
local coefmidfvn = results[1,1]
local semid = results[2,1]
local lowbmidfvn = results[5,1]
local upbmidfvn = results[6,1]
local coefmidfvncons = results[1,2]
local lowbmidfvncons = results[5,2]
local upbmidfvncons = results[6,2]
local nfvn = e(N)
****************
*90
logit policyadopted pred90_sw_lor if abs(pred90_sw - pred50_sw)>.1
*Store results as local vars
mat def results = r(table)
local coefhighfvn = results[1,1]
local sehigh = results[2,1]
local lowbhighfvn = results[5,1]
local upbhighfvn = results[6,1]
local coefhighfvncons = results[1,2]
local lowbhighfvncons = results[5,2]
local upbhighfvncons = results[6,2]





************************************************
************************************************
*Visualize Preference Gap Distribution in the Data Gilens' analysis creates
************************************************
graph set window fontface "Times New Roman"

******************************************
*FIGURE 1
******************************************
*90th income percentiles
******
*Change proportion variables to percent
gen pred90_sw_perc = 100*pred90_sw
gen pred10_sw_perc = 100*pred10_sw




*set number of bins 
local bins=40 
twoway ///
(histogram pred90_sw_perc if (pred10_sw_perc - pred90_sw_perc)>10, freq fcolor(navy%80) lcolor(black) lwidth(vvthin) bin(`bins')) ///
(histogram pred90_sw_perc if (pred90_sw_perc - pred10_sw_perc)>10, freq fcolor(orange%60) lcolor(black) lwidth(vvthin) bin(`bins')) ///
(pcarrowi 21 89 1.1 94, lcolor(black) mcolor(black)) ///
(pcarrowi 21 84.4 17 81, lcolor(black) mcolor(black)) ///
(pcarrowi 14 15 1.1 8, lcolor(black) mcolor(black)) ///
(pcarrowi 14 20 10.2 22, lcolor(black) mcolor(black)), ///
ylabel(, nogrid notick) xlabel(10(10)90, nogrid) yscale(lwidth(none)) /// 
xscale(lwidth(none) range(5 95)) ///
ytitle("Frequency", size(medlarge) axis(1)) ylabel(0(5)25, nogrid notick axis(1)) ///
yscale(lwidth(none) r(0 25)) ///
xtitle("Percent Favoring Policy Change (90th Income Percentile)", size(medlarge)) ///
plotregion(lcolor(black) lwidth(medium)) graphregion(fcolor(white) lcolor(white)) legend(off) ///
ysize(3) xsize(4) ///
text(24 87 "High income" "{it:more} supportive" "than low income" "(90th - 10th)>10%") ///
text(17 18 "High income" "{it:less} supportive" "than low income" "(10th - 90th)>10%") 

`FIG'graph export Figures/histogram90.eps, replace
`FIG'graph export Figures/histogram90.pdf, as(pdf) replace

******
*10th income percentiles
******
*set number of bins
local bins=40 
twoway ///
(histogram pred10_sw_perc if (pred10_sw_perc - pred90_sw_perc)>10, freq fcolor(navy%80) lcolor(black) lwidth(vvthin) bin(`bins')) ///
(histogram pred10_sw_perc if (pred90_sw_perc - pred10_sw_perc)>10, freq fcolor(orange%60) lcolor(black) lwidth(vvthin) bin(`bins')) ///
(pcarrowi 16.2 84 13.1 86.3, lcolor(black) mcolor(black)) ///
(pcarrowi 16.2 80 13.1 75.7, lcolor(black) mcolor(black)) ///
(pcarrowi 19.3 17 6.1 12.3, lcolor(black) mcolor(black)) ///
(pcarrowi 19.3 21 14.5 23.1, lcolor(black) mcolor(black)), ///
ylabel(, nogrid notick) xlabel(10(10)90, nogrid) yscale(lwidth(none)) /// 
xscale(lwidth(none) range(5 96)) ytitle("Frequency", size(medlarge)) /// 
xtitle("Percent Favoring Policy Change (10th Income Percentile)", size(medlarge)) ///
plotregion(lcolor(black) lwidth(medium)) graphregion(fcolor(white) lcolor(white)) legend(off) ///
ysize(3) xsize(4) ///
text(18.7 82 "High income" "{it:less} supportive" "than low income" "(10th - 90th)>10%") ///
text(21.7 19.8 "High income" "{it:more} supportive" "than low income" "(90th - 10th)>10%") 

`FIG'graph export Figures/histogram10.eps, replace
`FIG'graph export Figures/histogram10.pdf, as(pdf) replace
************************************************
************************************************


*************************************************************************************************
*Figure 2: Add predicted values
*************************************************************************************************
************************************************
*10th (90 vs 10)
************************************************

logit policyadopted pred10_sw_lor  if  abs(pred90_sw - pred10_sw)>.1
mat define b = get(_b)
gen yhat10=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred10_sw_lor  if (pred90_sw - pred10_sw)>.1
mat define b = get(_b)
gen yhat10_greater=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred10_sw_lor  if (pred10_sw - pred90_sw)>.1
mat define b = get(_b)
gen yhat10_lesser=exp(b[1,2]+(b[1,1]*pred10_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred10_sw_lor))) if e(sample)

local bins=40 
twoway ///
(histogram pred10_sw_perc if (pred10_sw_perc - pred90_sw_perc)>10, freq fcolor(navy%80) ///
lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1)) /// 
(histogram pred10_sw_perc if (pred90_sw_perc - pred10_sw_perc)>10, /// 
freq fcolor(orange%60) lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1) ///
ytitle("Frequency", size(medlarge) axis(1)) ylabel(, nogrid notick axis(1))) ///
(scatter yhat10 pred10_sw_perc, yaxis(2) msize(medsmall) mcolor(black) mlcolor(none) ///
ytitle("Predicted Probability of Policy Change", ///
size(mlarge) axis(2)) ylabel(0(.2)1, nogrid notick axis(2))) ///
(scatter yhat10_greater pred10_sw_perc, yaxis(2) msize(medsmall) mcolor(orange%80) mlcolor(none)) ///
(scatter yhat10_lesser pred10_sw_perc,  yaxis(2) msize(medsmall) mcolor(navy%90) mlcolor(none)),  ///
xtitle("Percent Favoring Policy Change (10th Income Percentile)", size(medlarge)) /// 
graphregion(fcolor(white) lcolor(white)) ///
plotregion(lcolor(black) lwidth(medium)) ///
xlabel(10[10]90, nogrid) xscale(lwidth(none) r(5 95)) ///
yscale(lwidth(none) r(0 1)) legend(off) ///
ysize(3) xsize(4) 

`FIG'graph export Figures/predicted10.eps, replace
`FIG'graph export Figures/predicted10.pdf, as(pdf) replace

************************************************
*90th (90 vs 10)
************************************************
logit policyadopted pred90_sw_lor  if  abs(pred90_sw - pred10_sw)>.1
mat define b = get(_b)
gen yhat90=exp(b[1,2]+(b[1,1]*pred90_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred90_sw_lor))) if e(sample)

logit policyadopted pred90_sw_lor  if (pred90_sw - pred10_sw)>.1
mat define b = get(_b)
gen yhat90_greater=exp(b[1,2]+(b[1,1]*pred90_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred90_sw_lor))) if e(sample)
mat drop b

logit policyadopted pred90_sw_lor  if (pred10_sw - pred90_sw)>.1
mat define b = get(_b)
gen yhat90_lesser=exp(b[1,2]+(b[1,1]*pred90_sw_lor))/(1+exp(b[1,2]+(b[1,1]*pred90_sw_lor))) if e(sample)
mat drop b

local bins=40 
twoway ///
(histogram pred90_sw_perc if (pred10_sw_perc - pred90_sw_perc)>10, freq fcolor(navy%80) ///
lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1)) ///
(histogram pred90_sw_perc if (pred90_sw_perc - pred10_sw_perc)>10, /// 
freq fcolor(orange%60) lwidth(vvthin) lcolor(black) lwidth(vvthin) bin(`bins') yaxis(1) ///
ytitle("Frequency", size(mlarge) axis(1)) ylabel(0(5)25, nogrid notick axis(1)) ///
yscale(lwidth(none) r(0 25))) ///
(scatter yhat90 pred90_sw_perc, yaxis(2) msize(medsmall) mcolor(black) mlcolor(none) ///
ytitle("Predicted Probability of Policy Change", ///
size(mlarge) axis(2)) ylabel(0(.2)1, nogrid notick axis(2))) ///
(scatter yhat90_greater pred90_sw_perc, yaxis(2) msize(medsmall) mcolor(orange%80) mlcolor(none)) ///
(scatter yhat90_lesser pred90_sw_perc, yaxis(2) msize(medsmall) mcolor(navy%90) mlcolor(none)), ///
xtitle("Percent Favoring Policy Change (90th Income Percentile)", ///
size(mlarge)) graphregion(fcolor(white) lcolor(white)) ///
plotregion(lcolor(black) lwidth(medium)) ///
xlabel(10[10]90, nogrid) xscale(lwidth(none) r(5 95)) ///
yscale(lwidth(none) r(0 1)) legend(off) ///
ysize(3) xsize(4) 

`FIG'graph export Figures/predicted90.eps, replace
`FIG'graph export Figures/predicted90.pdf, as(pdf) replace


******************************************
******************************************
*Figure 5
*Demonstrate responsiveness (beta) vs inercept (alpha)
******************************************
******************************************
*Identify specific policies
list YEAR QuestionText pred90_sw_perc pred10_sw_perc if pred90_sw_perc<24.02 & (pred90_sw - pred10_sw)>.1
list YEAR QuestionText pred90_sw_perc pred10_sw_perc if pred90_sw_perc>24.1 & pred90_sw_perc<24.2 & pred90_sw_perc!=. & (pred10_sw - pred90_sw)>.1
list YEAR QuestionText pred90_sw_perc pred10_sw_perc if pred90_sw_perc>17 & pred90_sw_perc<17.5 & pred90_sw_perc!=. & (pred10_sw - pred90_sw)>.1

list YEAR QuestionText pred90_sw_perc pred10_sw_perc if pred90_sw_perc>80 & pred90_sw_perc!=. & (pred10_sw - pred90_sw)>.1
list YEAR QuestionText pred90_sw_perc pred10_sw_perc if pred90_sw_perc>83.5 & pred90_sw_perc<83.7 & (pred90_sw - pred10_sw)>.1
list YEAR ID8102 QuestionText pred90_sw_perc pred10_sw_perc if pred90_sw_perc>81 & pred90_sw_perc<81.5 & (pred90_sw - pred10_sw)>.1

twoway scatter yhat90 pred90_sw_perc, msize(medsmall) mcolor(black) mlcolor(none) ///
|| scatter yhat90_greater pred90_sw_perc, msize(medsmall) mcolor(orange%80) mlcolor(none) ///
|| (scatter yhat90_lesser pred90_sw_perc, msize(medsmall) mcolor(navy%80) mlcolor(none)) ///
(pcarrowi .5 16.95851 .31 16.95851, lcolor(black) mcolor(black) lwidth(thick)) /// Cut medicare
(pcarrowi .1 20 .13 17.19092, lcolor(black) mcolor(black) lwidth(thick)) /// Increase defense budget
(pcarrowi .8 48 .34 24.01932, lcolor(black) mcolor(black)) /// Supply weapons to China
(pcarrowi .58 48 .16 24.14261, lcolor(black) mcolor(black)) /// Build fence (1993)
(pcarrowi .1 68 .385 81.4925, lcolor(black) mcolor(black) lwidth(thick)) /// state/federal college support
(pcarrowi .2 87 .48 81.38987, lcolor(black) mcolor(black) lwidth(thick)) /// eliminate marriage penalty
(pcarrowi .84 93 .49 83.62995, lcolor(black) mcolor(black)) /// line-item veto
(pcarrowi .48 98 .41 83.98826, lcolor(black) mcolor(black)), /// budget surplus medicare
graphregion(margin(r+30)) ///
ylabel(0(.2)1, nogrid notick) ylabel(, nogrid notick) /// 
xtitle("Percent Favoring Policy Change (90th Income Percentile)", size(medlarge)) ///
ytitle("Predicted Probability of Policy Change", size(medlarge)) ///
plotregion(lcolor(black) lwidth(medium)) xlabel(10[10]90, nogrid) xscale(lwidth(none) r(5 95)) yscale(lwidth(none) r(0 1)) legend(off) ///
ysize(3) xsize(5)  ///
text(.6 18 "Cut Medicare spending" "by $700 million" "({it:90th: 17% support}" "{it:10th: 6% support})", box bcolor(white)) ///
text(.06 34 "Increase defense budget" "({it:90th: 17% support}" "{it:10th: 30% support})") ///
text(.9 55 "Supply military weapons" "to Communist China" "({it:90th: 24% support}" "{it:10th: 12% support})") ///
text(.6 60 "Build a fence along entire" "border b/t U.S. & Mexico (1993)" "({it:90th: 24% support}" "{it:10th: 41% support})") ///
text(.08 80 "State or federal assistance for students with the" "ability but not enough money to attend college" "({it:90th: 81% support, 10th: 94% support})", box bcolor(white))  ///
text(.25 104 "Eliminating the so-called marriage" "penalty, which taxes a married couple" "({it:90th: 81% support, 10th: 52% support})", box bcolor(white))  ///
text(.48 110 "$600 billion of" "budget surplus over" "next 15 years to Medicare" "({it:90th: 84% support}" "{it:10th: 94% support})", box bcolor(white)) ///
text(.9 102 "Giving the President a line-item veto" "({it:90th: 84% support, 10th: 71% support})", box bcolor(white))

`FIG'graph export Figures/implications.eps, replace
`FIG'graph export Figures/implications.pdf, as(pdf) replace  
******************************************
******************************************
******************************************
******************************************

**********************************************************************************
**********************************************************************************
*export stored results to .txt file to insert in final manuscript (.tex file)
*Necessary packages
*ssc install `TEX'texresults2
**********************************************************************************
**********************************************************************************
*50 (50th vs 90th)
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coefmidfvn) result(`coefmidfvn') replace //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(semid) result(`semid') append //se
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowbmidfvn) result(`lowbmidfvn') append //lb 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upbmidfvn) result(`upbmidfvn') append //ub
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coefmidfvncons) result(`coefmidfvncons') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowbmidfvncons) result(`lowbmidfvncons') append //lb 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upbmidfvncons) result(`upbmidfvncons') append //ub
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(nfvn) result(`nfvn') round(0) append
*90 (50th vs 90th)
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coefhighfvn) result(`coefhighfvn') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(sehigh) result(`sehigh') append //se
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowbhighfvn) result(`lowbhighfvn') append //lb
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upbhighfvn) result(`upbhighfvn') append //ub
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coefhighfvncons) result(`coefhighfvncons') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowbhighfvncons) result(`lowbhighfvncons') append //lb 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upbhighfvncons) result(`upbhighfvncons') append //ub
*10 (10th vs 90)
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coeflowtvn) result(`coeflowtvn') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(selow) result(`selow') append //se
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowblowtvn) result(`lowblowtvn') append //lb
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upblowtvn) result(`upblowtvn') append //ub 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coeflowtvncons) result(`coeflowtvncons') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowblowtvncons) result(`lowblowtvncons') append //lb 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upblowtvncons) result(`upblowtvncons') append //ub
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(ntvn) result(`ntvn') round(0) append
*90 (10th vs 90th)
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coefhightvn) result(`coefhightvn') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(sehightvn) result(`sehightvn') append //se
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowbhightvn) result(`lowbhightvn') append //lb
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upbhightvn) result(`upbhightvn') append //ub
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(coefhightvncons) result(`coefhightvncons') append //coef
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(lowbhightvncons) result(`lowbhightvncons') append //lb 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(upbhightvncons) result(`upbhightvncons') append //ub
*Correlation: income group preferences
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(corlowhigh) result(`corlowhigh') round(2) append 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(cormidhigh) result(`cormidhigh') round(2) append 
`TEX'texresults2 using TxtFiles_NumericalResults\results.txt, texmacro(cormidpublic) result(`cormidpublic') round(2) append 
**********************************************************************************************
**********************************************************************************************






